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Abstract 

An  experimental  and  numerical  research  has  been  performed  in  order  to  study  the  flow  distribution  in  a  bipolar  plate  of  a  commercial  PEM 
fuel  cell.  Planar  laser  induced  fluorescence  (PLIF)  trace  tracking  has  been  applied  to  visualize  the  flow  pattern  and  to  measure  the  velocity 
in  the  plate  channels.  Simultaneously,  the  problem  has  been  studied  numerically,  simulating  the  flow  under  similar  operational  conditions 
as  those  fixed  in  the  experiments.  Results  obtained  reveal  a  defective  design  of  the  bipolar  plate.  Based  on  the  experimental  visualization 
and  on  the  numerical  simulations  it  is  concluded  that  the  flow  preferentially  moves  through  the  lateral  channels,  resulting  in  an  inappropriate 
distribution  on  the  electrode  surfaces.  Velocity  measurements  also  confirm  the  above  statements,  showing  high  values  at  the  lateral  channels, 
while  the  flow  is  nearly  stagnant  in  the  central  region.  With  this  non-homogeneous  flow  distribution  at  the  bipolar  plate,  a  low  performance 
of  the  fuel  cell  energy  conversion  could  be  expected. 

©  2005  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

The  energy  demand  worldwide  is  growing  at  an  alarming 
rate.  This  demand  is  responded  by  an  increase  in  the  com¬ 
bustion  of  fossil  fuels,  with  the  entailed  problems  of  pollutant 
emissions,  greenhouse  effect  and  acid  rain.  Besides,  the  nat¬ 
ural  reserves  of  fossil  fuels  are  diminishing  and  a  substantial 
increase  in  their  price  can  be  expected  in  the  foreseeable  fu¬ 
ture.  In  view  of  this  situation,  hopes  have  been  deposited 
in  fuel  cells  as  a  key  solution  for  the  21st  century  energy 
problems,  enabling  clean  and  efficient  production  of  heat  and 
power  from  a  diversity  of  primary  sources  [1]. 

Fuel  cells  are  devices  that  generate  electricity  by  a  chem¬ 
ical  reaction.  Despite  the  general  belief  that  they  represent 
a  new  technology,  their  basic  working  principles  have  been 
known  for  centuries.  Strictly  speaking,  they  were  first  demon¬ 
strated  in  1 843  with  the  experiments  of  Sir  William  R.  Grove. 
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Fuel  cells  have  two  electrodes,  anode  (negative)  and  cathode 
(positive),  where  the  reaction  that  produces  electricity  takes 
place.  Between  them,  there  is  an  electrolyte,  which  allows  the 
transport  of  electrically  charged  particles  from  one  electrode 
to  the  other,  and  a  catalyst,  which  enables  and  speeds  up  the 
reaction.  Depending  on  the  electrolyte,  fuel  cells  are  usually 
divided  into  five  different  types,  namely:  alkali,  molten  car¬ 
bonate,  phosphoric  acid,  proton  exchange  membrane  (PEM) 
and  solid  oxides.  In  general,  hydrogen  is  normally  used  as 
the  basic  fuel,  which  reacts  electrochemically  with  oxygen 
producing  electricity  and  heat  and,  as  a  byproduct,  water. 

Proton  exchange  membrane  fuel  cells  work  with  a  poly¬ 
meric  electrolyte  formed  into  a  thin,  very  light  and  permeable 
sheet,  at  low  temperature  (around  80  °C).  To  speed  up  the 
chemical  reaction,  a  platinum  catalyst  is  used  in  both  sides 
of  the  membrane.  Electrons  are  removed  from  the  hydrogen 
atoms  in  the  anode  and  the  resulting  protons  diffuse  through 
the  permeable  membrane,  migrating  towards  the  cathode. 
The  electrons  pass  from  the  anode  to  the  cathode  through  an 
external  circuit  and  provide  electric  power.  At  the  cathode, 


F.  Barreras  et  al.  /  Journal  of  Power  Sources  144  (2005)  54-66 


55 


Nomenclature 

d  inlet  duct  diameter  (m) 

E  error  in  the  numerical  simulation 

n  iteration  number  in  the  numerical  study 

Q  volumetric  flow  rate  (ml  min- 1  or  1  min- 1 ) 

Re  Reynolds  number 

v  velocity  of  the  flow  (ms-1  or  cm s  - 1 ) 

V  volume  (m3) 

X  fraction  of  liquids  in  the  mixture  (%) 

Greek  symbols 

P  criteria  of  convergence  in  the  numerical  simu¬ 

lation 

0  solution  of  the  problem  in  the  numerical  study 

fi  absolute  viscosity  (Pa  s) 

p  density  (kg  m- 3 ) 

v  kinematics  viscosity  (m2  s-1) 

Subscripts  and  superscripts 
glyc.  glycerin 

H2  hydrogen 

liq.  liquids 

mixt.  mixture 

wat.  water 


electrons,  protons  and  oxygen  combine  to  form  water.  PEM 
fuel  cells  have  been  widely  used  since  the  early  days  of  the 
space  programs,  as  well  as  in  submarine  vessels.  The  strict 
environmental  regulations  on  air  quality  have  also  motivated 
the  introduction  of  this  technology  in  the  automotive  indus¬ 
try.  In  the  same  way,  PEM  fuel  cells  have  also  been  used  in 
stationary  power  applications.  Today,  the  broad  spectrum  of 
applications  for  PEM  fuel  cells  includes  their  use  in  portable 
devices  such  as  cellular  telephones  or  laptop  computers. 

The  efficiency  in  energy  conversion  achieved  in  PEM  fuel 
cells  is  higher  than  that  in  power  plants  or  internal  combustion 
engines.  This  efficiency  can  be  reached  as  the  result  of  a  set  of 
complex  physical  and  chemical  processes  occurring  simulta¬ 
neously,  which  are  strongly  dependent  on  the  fuel  and  oxygen 
fluid  dynamics  inside  the  fuel  cell.  It  is  easy  to  understand  that 
optimization  of  PEM  fuel  cells  performance  requires  a  deep 
comprehension  of  the  current  density  behavior  as  a  function 
of  the  operational  conditions.  In  this  sense,  several  studies 
have  been  performed  applying  numerical  techniques  [2-4]. 
From  the  experimental  viewpoint,  some  methods  used  have 
been  aimed  to  study  the  water  production  and  its  influence 
on  the  mass  transport  inside  the  fuel  cell  [5-9] . 

A  singular  component  of  PEM  fuel  cells  is  the  bipolar 
plate,  as  the  first  stage  of  the  flow  distribution  system.  In 
general,  the  functions  of  this  element  are  to  guide  the  flow 
of  reactant  gases,  and,  in  the  case  of  fuel  cell  stacks,  to  pro¬ 
vide  electron  conduction  paths  from  one  cell  to  another.  So, 


in  order  to  ensure  a  correct  operation  of  the  fuel  cell  it  is 
necessary  both  a  suitable  material  and  the  proper  topology  to 
optimize  the  gas  flow  distribution.  Common  materials  used 
in  bipolar  plates  are  graphite,  some  metals  (stainless  steel, 
titanium,  etc.),  or  composite  materials  [10],  which  should 
be  lightweight,  and  easy  to  manufacture.  They  must  also  be 
impermeable  to  gases,  and  have  high  electrical  and  thermal 
conductivities.  A  low  performance  of  the  fuel  cell  energy 
conversion,  as  well  as  a  waste  of  the  very  expensive  plat¬ 
inum  catalyst  could  be  expected  if  the  flow  distribution  in  the 
bipolar  plate  is  non-homogeneous. 

There  are  different  plate  structures  that  produce  different 
flow  fields  as,  for  example,  squared  spot,  interdigitated,  ser¬ 
pentine,  spirals,  cascade  and  series-parallel,  with  geometries 
of  various  complexity  levels.  Alternatively,  flow  through 
porous  carbon  or  perforated  stainless  steel  plates  have  also 
been  used.  A  comparative  analysis  of  the  advantages  and 
disadvantages  of  some  of  the  above  topologies  can  be  found 
in  the  reviews  by  Carrette  et  al.  [11]  or  Costamagna  and 
Srinivasan  [12].  Hontanon  et  al.  [13]  studied  the  gas  flow 
distribution  system  of  a  PEM  fuel  cell  using  3D  numerical 
simulations.  Their  results  indicated  that  porous  materials 
yielded  better  flow  distributions  compared  to  grooved  plates 
in  terms  of  reactant  gas  utilization.  Recently,  Dohle  et  al. 
[14]  studied  the  interaction  between  the  diffusion  layer 
and  the  flow  field  exiting  a  meander  channel  bipolar  plate 
in  a  direct  methanol  fuel  cell.  Based  also  on  numerical 
simulations,  they  concluded  that  even  a  meander  struc¬ 
ture  distributes  the  reactants  non-homogeneously  on  the 
electrodes. 

In  this  study,  an  experimental  and  numerical  research  has 
been  performed  in  order  to  study  the  flow  distribution  in  a 
commercial  parallel  channel  bipolar  plate.  To  this  end,  flow 
visualization  using  laser-induced  fluorescence,  as  well  as 
measurements  of  the  velocity  field  by  dye  trace  tracking  have 
been  applied.  On  the  other  hand,  a  2D  numerical  simulation 
of  the  flow  distribution  based  on  the  Navier-Stokes  equations 
has  also  been  performed.  Results  obtained  from  both  exper¬ 
imental  and  numerical  studies  have  been  compared.  Both  of 
them  reveal  a  non-homogeneous  flow  distribution  across  the 
bipolar  plate,  which  will,  probably  produce  a  limited  per¬ 
formance  of  the  fuel  cell  energy  conversion.  The  excellent 
agreement  between  the  experimental  results  and  the  numer¬ 
ical  predictions  confirms  the  validity  of  the  numerical  code 
to  study  design  variations  without  the  need  of  actually  fabri¬ 
cating  the  plates. 

2.  The  problem  under  study 

The  experiments  have  been  performed  using  a  commer¬ 
cial  bipolar  plate  with  a  rough  surface  area  of  50  cm2.  It  is 
formed  by  16  central  channels  1  mm  deep  and  3  mm  wide, 
separated  by  ribs  with  a  1  mm  thickness.  Two  channels  in  the 
upper  and  lower  part  of  the  plate  and  two  lateral  channels,  all 
of  them  with  a  width  of  2  mm,  surround  the  whole  flow  area. 
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The  internal  diameter  of  both  inlet  and  outlet  ducts  is  2  mm. 
This  type  of  bipolar  plates  is  normally  used  to  distribute  the 
oxygen  (or  air)  and  hydrogen  over  the  cathode  and  anode, 
respectively.  However,  this  study  has  been  focused  on  the  vi¬ 
sualization  of  the  hydrogen  flow  pattern  distribution.  In  this 
case,  a  single  fuel  cell  with  a  volumetric  flow  of  1  lmin-1 
(1.67  x  1(T5  m3  s_1)  at  80  °C  has  been  considered  [15].  Un¬ 
der  this  flow  condition,  the  single  fuel  cell  should  produce  a 
power  close  to  50  W,  corresponding  to  a  current  density  of 
400  mA  cm-2  and  a  voltage  of  2.5  V. 

Gas  flow  visualization  has  been  widely  used  in  fluid  me¬ 
chanics  but  its  application  involves  some  complexities  diffi¬ 
cult  to  circumvent  in  some  experimental  setups.  For  exam¬ 
ple,  ultraviolet  laser  radiation  is  needed  to  induce  the  flu¬ 
orescence  of  the  typical  vapor  tracers  used  to  visualize  the 
flows,  such  as  acetone  or  acetaldehyde.  In  this  case,  a  quartz 
window  is  needed  to  seal  the  bipolar  plate  in  the  front-face 
where  the  images  are  to  be  recorded.  Additionally,  these  or¬ 
ganic  vapors  cannot  be  used  in  conjunction  with  many  plastic 
materials  that  are  quickly  degraded.  For  these  reasons,  a  di¬ 
mensional  analysis  has  been  performed  to  replace  the  gas  by 
a  liquid,  ensuring  the  dynamic  similarity  between  the  real 
problem  and  the  model.  The  main  physical  properties  of  hy¬ 
drogen  at  80  °C  considered  in  the  analysis  [16,17]  are:  abso¬ 
lute  viscosity  (/x)  of  1.05  x  10-5  Pas  and  density  ( p )  equal 
to  0.0674  kg  m-3,  resulting  in  a  kinematic  viscosity  (v)  of 
1.558  x  10-4  m2  s-1. 

A  gas  entrance  velocity  of  5.3  m  s_1  results  from  a  simple 
calculation  for  the  hydrogen  flow  ( vu2 )  using  the  volumetric 
flow  rate  considered  and  the  area  of  the  inlet  duct.  The  result¬ 
ing  Reynolds  number  (Re  =  vu 2d/v)  based  on  the  inlet  duct 
diameter  (d)  is  68.  In  order  to  perform  the  experiments  under 
dynamic  similarity  when  a  liquid  is  used,  new  operational 
conditions  have  to  be  considered  to  maintain  the  Reynolds 
number  unchanged,  as  this  dimensionless  group  is  the  one 
that  essentially  determines  the  flow  characteristics.  Two  con¬ 
ditions  are  imposed  to  the  selected  liquid:  it  has  to  be  transpar¬ 
ent,  and  the  fluorescent  dye  that  will  be  used  has  to  be  easily 
soluble  in  it.  For  simplicity,  only  three  candidates  have  been 
considered,  water,  glycerin  and  methanol.  The  main  physi¬ 
cal  properties  for  the  three  fluids  at  20  °C  [16]  are  shown  in 
Table  1 .  If  flow  dynamic  similarity  is  imposed  to  the  selected 


liquid,  Eq.  (1)  is  used  to  calculate  the  liquid  velocity,  as  well 
as  the  volumetric  flow  rate 

Uia  d 

Re  mod  =  68  =  (1) 

Uiq. 

Results  obtained  applying  Eq.  (1)  to  the  three  selected 
liquids  are  displayed  in  Table  2.  It  is  evident  that  extremely 
low  flow  velocities  have  to  be  managed  if  water  or  methanol  is 
used  in  the  experiments.  Under  these  experimental  conditions 
the  liquid  flow  will  be  largely  perturbed  by  the  dye  injection. 
On  the  other  hand,  a  quite  high  velocity  flow  has  to  be  supplied 
if  glycerin  is  used,  due  to  its  high  viscosity  value.  In  this 
case,  a  huge  pressure  is  needed  to  overcome  the  pressure 
drop  inside  the  bipolar  plate,  making  the  watertight  seal  of 
the  model  a  difficult  task. 

To  solve  all  of  these  problems  a  mixture  of  glycerin  and 
water  has  been  selected.  For  example,  if  a  volumetric  flow 
rate  of  230  ml  min-1  is  fixed,  which  corresponds  to  a  velocity 
of  1.22  ms-1,  a  mixture  kinematic  viscosity  value  of  24.1 
cSt(24.1  x  10  6m2s  l)  is  needed.  The  mixture  composition 
can  be  calculated  using  the  methodology  described  in  [16]  for 
aqueous  systems  applying  the  following  set  of  equations: 

XWat  ^glyc. 

/WMglyc. 

Unixt.  —  yJ') 

Anixt. 

Anixt.  —  2fwat.  Pwat.  T  ^glyc.Pglyc.  (3) 


where  X  represents  the  volumetric  fraction  of  both  liquids 
defined  by 


kjiq. 


Solving  the  equations,  a  kinematic  viscosity  of 
24.1  x  10  6  m2  s  1  requires  a  mixture  composition  of  49.6% 
of  glycerin  and  5 1 .4%  of  water. 


3.  Experimental  facilities 

Having  in  mind  the  above-calculated  values  for  the  ex¬ 
perimental  conditions,  the  experimental  facility  displayed  in 
Fig.  1  has  been  used  in  the  experiments.  It  includes  the  bipolar 


Table  1 


Physical  properties  of  the  liquids  considered  in  the  study  at  20  °C 


Hydrogen 

Water 

Glycerin 

Methanol 

Density,  p  (kgm-3) 

Absolute  viscosity,  /x  (Pa  s) 
Kinematic  viscosity,  v  (m2  s-1) 

0.0674 

1.05  x  10"5 

1.558  x  10“4 

998 

1.003  x  10“3 

1.005  x  10“6 

1.236 

0.7601 

6.15  x  10-4 

788.5 

5.99  x  10-4 
7.6  x  10“7 

Table  2 

Velocity  and  volumetric  flow  rate  calculation  for  the  three  liquid  fluids  considered 

Hydrogen 

Water 

Glycerin 

Methanol 

Velocity,  v 

Volumetric  flow  rate,  Q 

5.3ms_1 

1.01  min-1 

3.42  cms-1 

6.45  ml  min-1 

20.91  ms-1 

3.941  min-1 

2.58  cms-1 
4.86  ml  min-1 
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dye  injection  point 


Fig.  1.  Experimental  setup. 


plate,  a  101  tank,  an  immersion  pump  capable  of  supplying 
a  maximum  flow  rate  of  1001  min-1,  two  manometers  and 
a  flowmeter  with  a  maximum  reading  of  1.7 1  min-1.  A  by¬ 
pass  has  been  included  before  the  bipolar  plate  assembly  to 
avoid  large  pressure  drops  and,  hence,  liquid  overheating. 
This  is  very  important  due  to  the  exponential  dependence 
of  the  liquid  mixture  viscosity  on  temperature,  which  might 
cause  large  fluctuations  on  the  flow  conditions.  The  dye  in¬ 
jection  point  is  placed  80  mm  (40  diameters)  upstream  of  the 
inlet  point,  and  consists  of  a  1  mm  capillary  welded  at  an 
angle  of  45°  with  respect  to  the  axis  of  the  inlet  duct.  A  pho¬ 
tograph  of  the  injection  system  can  be  seen  in  Fig.  2.  This 
figure  also  shows  the  transparent  plastic  cover  that  seals  the 
bipolar  plate  and  constitutes  the  front-face  of  the  experimen¬ 
tal  model.  It  is  to  be  noted  that  the  inlet  and  outlet  ports  are 
tilted  with  respect  to  the  plate  edges. 

Two  different  flow  rates  have  been  tested,  230  and 
350  ml  min-1.  Before  acquiring  the  images,  the  liquid  was 
allowed  to  circulate  during  5  min  to  stabilize  the  flow  in  sta¬ 
tionary  conditions.  After  this  time,  dye  was  injected  in  pulses, 
in  synchronization  with  the  laser  pulses  and  the  image  acqui¬ 
sition. 

To  visualize  the  flow  pattern,  laser  induced  fluorescence 
(LIF)  has  been  applied,  using  Sulforhodamine  B  (Kiton  Red) 
dissolved  in  the  water/glycerin  mixture  as  the  luminescent 
tracer.  To  excite  the  dye,  a  double  cavity  Quantel  YG78 1C-10 


Fig.  2.  Detail  of  the  dye  injection  system. 


pulsed  Nd:  YAG  laser  has  been  used,  doubling  the  frequency 
of  its  emission  to  obtain  100  mJ  pulses  at  532  nm,  with  a  pulse 
duration  of  6  ns.  The  absorption  spectrum  of  Sulforhodamine 
B  has  a  maximum  at  556  nm;  hence  excitation  at  532  nm  is 
very  efficient  [18].  The  fluorescence  peak  of  this  dye  is  lo¬ 
cated  at  620  nm,  which  allows  for  an  efficient  discrimination 
between  the  emitted  signal  and  the  excitation  laser  light  beam. 
To  further  decrease  the  background  light,  a  Schott  OG  550 
filter  has  been  placed  in  front  of  the  camera  lens,  blocking  any 
residual  532  nm  wavelength  from  the  incident  laser  beam. 

To  image  the  fluorescence  emission,  a  Princeton  Instru¬ 
ments  slow  scan  CCD  camera  has  been  used,  with  a  50  mm 
/#1.2  Nikon  lens,  placed  perpendicular  to  the  bipolar  plate. 
As  depicted  in  Fig.  3,  a  spherical  lens  with  a  focal  distance  of 
—25  mm  has  been  placed  at  1  m  of  the  bipolar  plate  to  expand 
the  laser  beam  obtaining  a  slightly  divergent  cylinder  of  light, 
which  illuminates  the  plate  at  an  angle  of  60°  with  respect 
to  the  axis  of  the  camera.  This  arrangement  does  not  intro¬ 
duce  any  perspective  distortion  in  the  images  and  eliminates 
any  possible  direct  reflection  of  the  laser  light  on  the  CCD 
camera.  Data  sets  have  been  recorded  for  a  field  of  view  of 
70  mm  x  94  mm  with  a  spatial  resolution  of  245  |jim  pixel- 1 . 

Velocity  field  in  the  channels  has  been  measured  by  track¬ 
ing  the  evolution  of  the  dyed  liquid  using  the  same  optical 


spherical  lens/=-25  mm  @  532  nm  laser  beam 


Fig.  3.  Optical  setup  for  visualization  and  PIV  experiments. 
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arrangement  and  laser  source  previously  described  for  the  vi¬ 
sualization  experiments.  An  interlined  Hamamatsu  ORCA- 
ER  CCD  camera  with  a  resolution  of  1024  x  1024  pixels  has 
been  employed  to  record  the  image  sequences.  The  time  in¬ 
terval  between  successive  images  was  set  to  0.5  s  and  the 
correlation  between  them  has  been  analyzed  to  track  the  dis¬ 
placement  of  the  dyed  liquid  front. 


4.  Numerical  simulation 

In  order  to  numerically  study  the  present  configuration, 
a  computer  code  has  been  adapted  to  simulate  the  flow  de¬ 
scribed  in  the  previous  section.  As  numerical  simulations  are 
relatively  cheap,  this  code  will  be  a  tool  to  study  possible 
alternative  geometrical  configurations  in  the  future. 

A  first  approach  to  simulate  this  flow  is  to  consider  a  2D 
geometry.  In  principle  this  seems  a  rather  limiting  simplifi¬ 
cation,  given  the  non-homogeneous  nature  of  the  flow  in  the 
coordinate  removed.  However,  the  most  important  issue  that 
influences  the  flow  pattern  is  the  pressure  distribution  near  the 
entrance.  This  pressure  distribution  can  be  adequately  repro¬ 
duced  with  a  2D  simulation,  so  the  flow  rate  in  each  channel  of 
the  plate  will  match  reasonably  well  the  experimental  results, 
as  it  will  be  seen  later.  Other  assumptions,  such  as  isothermal, 
incompressible  and  laminar  flow,  and  a  steady  final  solution, 
have  been  included  to  match  the  experimental  conditions. 

The  resulting  transport  equations  governing  the  flow  in  the 
channels  of  the  bipolar  plate  are  the  steady  version  of  the  well 
known  Navier-Stokes  equations  for  mass  and  momentum 
conservation: 

V(pu)  =  0  (5) 

V(pv  •  n)  =  —V p  +  V(/xV t>)  +  pg  (6) 

A  finite  volume  method  [19]  has  been  preferred  over  finite 
elements  to  discretize  the  coupled  Eqs.  (5)  and  (6),  given 
the  simple  geometrical  configuration  of  the  channels  in  the 
bipolar  plate.  Finite  differences  have  not  been  considered 
due  to  the  non-conservative  character  of  this  discretization 
method. 

To  improve  the  velocity-pressure  coupling  a  (regular 
Cartesian)  staggered  grid,  that  is,  shifted  half  a  cell  in  v 
and  y  directions  for  calculating  v  and  y  components  of  the 
velocity  field,  has  been  used.  This  coupling  is  numerically 
implemented  using  a  SIMPLE  [20]  algorithm.  Finally,  a 
gradient-conjugate  technique  has  been  used  to  iteratively 
solve  the  resulting  linear  equation  system. 

To  accommodate  the  2D  numerical  simulation  to  the  3D 
physical  problem  is  essential  to  preserve  the  Reynolds  num¬ 
ber  in  the  channel.  For  a  real  flow  rate,  this  imposes  a  2D  flow 
rate  given  through  the  following  relationship: 

03D  =  -d02D.  (7) 


So,  for  a  volumetric  flow  of  230  ml  min-1  and  a  Reynolds 
number  of  68,  the  corresponding  2D  inlet  velocity  is 
1 .22  m  s_  1 ,  which  equals  that  obtained  from  the  experimental 
analysis. 

In  the  real  configuration,  the  inlet  duct  forms  an  angle 
with  the  XY  plane,  which  is  not  reproducible  in  a  2D  approx¬ 
imation.  As  a  result,  the  eddy  that  appears  near  the  entrance 
might  be  shifted  with  respect  to  the  experimental  position. 
This  eddy  can  be  influential  in  the  general  flow  pattern  in  the 
channel,  for  this  reason,  it  is  very  important  to  reproduce  its 
actual  position  as  accurately  as  possible.  This  has  been  done 
by  introducing  the  fluid  at  5  mm  from  the  corner  of  the  plate, 
and  by  properly  choosing  the  incidence  angle  (70°  respect 
to  the  horizontal  axis).  However,  it  has  been  verified  that  the 
actual  position  of  the  eddy  is  not  too  sensitive  to  this  two 
parameters. 

Although  the  flow  is  steady,  the  numerical  scheme  requires 
an  initial  condition  to  begin  the  simulation.  A  constant  (in 
modulus)  small  velocity  has  been  considered  through  all  the 
channels.  The  numerical  value  has  been  chosen  to  preserve 
mass  conservation. 


5.  Error  analysis 

To  study  the  effect  of  the  discretization  on  the  numerical 
solution,  different  grid  sizes  (69  x  69, 138  x  138, 207  x  207) 
have  been  considered.  The  calculation  of  the  velocity  field  of 
the  flow  for  these  three  grids  is  depicted  in  Fig.  4.  As  can  be 
seen,  the  flow  behavior  obtained  is  quite  similar  even  when 
the  flow  is  simulated  using  relatively  coarse  grids.  It  can  also 
be  observed  that  meshes  138  x  138  and  207  x  207  show  very 
similar  results,  which  indicates  that  the  intermediate  mesh 
has  a  spatial  resolution  high  enough  to  simulate  with  a  good 
precision  the  flow  conditions  in  the  bipolar  plate.  The  results 
of  the  69  x  69  mesh  converged  after  50,000  iterations  with 
a  level  of  convergence  /3  less  than  10-6,  being  even  lower 
for  the  other  two  grids.  This  coefficient  /3  is  defined  in  the 
following  way:  assuming  that  the  solution  after  each  iteration 
behaves  as  a  Cauchy  series,  an  error  can  be  defined  associated 
to  each  iteration  as 

Eo  =  101  -  00 1  >  =  |02  -  01 1  >  •  •  •  >  En 

=  1071+1  —  0n  I  >  •  •  •  0.  (8) 

Extending  the  previous  idea  to  define  an  error  associated 
between  two  levels  of  iteration,  using  the  triangle  inequality, 

10712  _  0711  I  —  1 0722  _  0712-1  H - 1-  10711  _  0711-1  I 

<  I n2  -  n\\\ 02  -  01 1,  (9) 

it  is  obtained 
^  _  10712  _  0711 

712— Til  = 


«2  -  nil  102  -  01 1 


(10) 
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Fig.  4.  Calculation  of  velocity  values  of  the  flow  in  the  different  channels  for  the  three  grids. 
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This  is  a  non-dimensional  error  with  values  between  0  and 
1 .  A  criteria  or  level  of  convergence  /3  is  reached  when 

En2-ni<P\  V*  (11) 

The  influence  of  the  number  of  iterations  in  the  conver¬ 
gence  of  the  solution  has  also  been  analyzed.  Results  are 
presented  in  Fig.  5  for  the  intermediate  grid  (138  x  138), 


which  shows  that  the  velocity  in  the  penultimate  channel  takes 
longer  to  converge.  A  physical  explanation  to  this  issue  is  that 
the  flow  in  this  channel  is  mainly  affected  for  the  presence  of 
the  outlet  duct,  which  is  very  close,  causing  a  slower  setting 
of  a  proper  pressure  field  at  its  exit.  To  minimize  errors  due 
to  this  influence,  the  first  5  mm  of  the  exit  pipe  have  been 
included  in  the  simulation  domain.  A  detailed  description  of 
the  convergence  error  for  this  grid  size  is  depicted  in  Fig.  6, 
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Fig.  5.  Influence  of  the  number  of  iterations  in  the  convergence  of  the  solution  of  the  numerical  study. 


Fig.  6.  Influence  of  the  convergence  error  for  the  138  x  138  grid  as  a  function  of  the  number  of  iteration  steps. 
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(a) 


(b) 


Fig.  7.  Visualization  images  obtained  for  two  different  water  flows:  (a)  230  ml  min 


and  (b)  350  ml  min 


which  shows  the  error  En2-ni  in  the  channel  cross  section 
averaged  velocities  for  six  different  number  of  iterations.  It 
certainly  shows  a  good  numerical  behavior. 

6.  Results 

As  a  first  step  of  the  experimental  research,  some 
preliminary  tests  were  performed  using  water.  The  aim  of 
these  experiments  was  to  check  the  overall  efficiency  of  the 
experimental  set  up.  In  this  case,  a  thorough  control  of  the 
temperature  variation  of  the  fluid  at  the  reservoir  to  evaluate 
the  efficiency  of  the  by-pass  system,  as  well  as  the  control  of 
the  seeding  method  was  performed.  In  Fig.  7,  a  sample  of  the 


acquired  images  can  be  observed  for  two  different  volumetric 
flow  rates,  230  and  350  ml  min-1,  which  correspond  to 
Reynolds  numbers  of  2. 120  and  3.7 14,  respectively.  It  should 
be  noted  that  in  all  the  images,  the  fluid  inlet  is  located  in 
the  upper  right  corner,  and  liquid  flows  from  right  to  left. 

Once  the  experimental  system  had  been  tested,  experi¬ 
ments  with  the  glycerin-water  mixture  were  performed.  The 
kinematic  viscosity  of  the  final  mixture  was  experimentally 
measured  with  a  Rheotest  viscosimeter,  adjusting  in  situ  the 
desired  value  by  slightly  varying  the  fraction  of  liquids  in 
the  mixture.  An  example  of  the  recorded  images  is  displayed 
in  Fig.  8  for  the  same  volumetric  flow  rates  as  in  Fig.  7,  230 
and  350  ml  min-1 ,  corresponding  to  Reynolds  numbers  of  68 
and  154.  One  striking  feature  of  Figs.  7  and  8  is  the  similar- 


Fig.  8.  Visualization  images  for  two  different  glycerin  +  water  flows:  (a)  230  ml  min 


and  (b)  350  ml  min 
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ity  of  the  dye  pattern  for  water  and  the  mixture  inside  the 
bipolar  plate.  The  large  difference  in  Reynolds  number  for 
both  flows  indicates  that  the  flow  distribution  is  independent 
either  of  the  flow  conditions  or  the  fluid  physical  properties 
and  suggests  very  low  turbulence  levels.  Flow  visualization 
results  show  a  very  poor  efficiency  of  the  bipolar  plate  tested 
regarding  to  flow  distribution.  The  presence  in  the  images  of 
channels  practically  depleted  of  any  dye  trace  is  indicative 
of  absence  of  flow  circulation,  in  contrast  to  the  preferred 


circulation  along  the  lateral  channels.  The  flow  of  dyed  fluid 
through  the  lateral  (upper  and  lower)  channels,  as  well  as  the 
higher  dilution  of  dye  concentration  in  these  areas  clearly  in¬ 
dicates  a  substantially  higher  velocity  compared  to  that  of  the 
central  zone  of  the  bipolar  plate  where  the  dye  concentration 
is  homogeneously  distributed. 

Fig.  9  presents  a  sequence  of  the  temporal  evolution  of 
the  dye  front  progressing  across  the  plate.  The  temporal  in¬ 
terval  between  successive  images  is  0.5  s.  Some  details  have 


Fig.  9.  Temporal  sequence  showing  the  advance  of  the  dyed  fluid.  Images  are  separated  by  0.5  s.  Glycerin  +  water  flow  rate  is  230  ml  min 
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Fig.  10.  Close  up  detail  of  the  recirculation  bubble  at  the  entrance  of  the 
central  channels. 

to  be  outlined.  It  is  clearly  observed  that  a  recirculation  bub¬ 
ble  forms  in  the  entrance  of  each  channel,  in  particular  in 
the  central  zone  of  the  plate  (Fig.  10).  This  feature  substan¬ 
tially  decreases  the  free  inlet  area  for  these  channels,  and 
contributes  to  explain  the  reduced  liquid  entrance  in  them. 
The  lower  velocity  in  the  third  channel,  compared  to  the  one 
immediately  below  it  is  partly  due  to  the  angle  that  the  injec¬ 
tion  port  forms  with  the  plate.  As  a  matter  of  fact,  the  flow 
along  the  central  channels  is  so  slow  that  in  certain  cases  it  is 
possible  to  observe  liquid  entering  through  the  exit  (Fig.  11). 

Numerical  simulation  also  shows  similar  results,  and  with 
a  good  agreement  with  the  experiments.  Following  the  thor¬ 
ough  discussion  in  the  Section  5,  the  intermediate  mesh 
(138  x  138)  has  enough  precision  to  obtain  reliable  results 
in  the  numerical  study.  For  this  reason,  all  the  results  pre¬ 
sented  here  are  referred  to  this  mesh,  unless  otherwise  is 
pointed  out.  In  Fig.  12,  comparisons  with  the  experiments 
are  shown  for  a  level  of  convergence  of  order  10-20  for 
^50,000—20,000  ?  which  seems  very  accurate,  taken  into  account 
that  a  2D  approximation  has  being  used.  In  order  to  ease 
the  comparison,  for  the  numerical  simulation  a  flow  rate  of 


Fig.  11.  Close  up  detail  of  the  plate  outlet. 


230  ml  min-1  of  glycerin-water  mixture  has  been  consid¬ 
ered,  which  corresponds  to  one  of  those  used  in  the  exper¬ 
iments.  The  agreement  with  the  experimental  data  is  quite 
noticeable  and  certainly  corroborates  the  poor  performance 
of  the  design  for  the  bipolar  plate  studied  in  the  present  work. 
Again,  it  has  been  verified  that  the  flow  circulates  with  a 
higher  velocity  through  the  upper,  lower  and  lateral  channels. 
Their  corresponding  pressure  and  velocity  maps  are  shown 
in  Figs.  13  and  14,  respectively.  Numbers  at  the  two  sides 
of  Fig.  13  correspond  to  the  pressure  values  at  the  entrance 
and  exit  of  each  channel,  respectively.  It  may  be  seen  that 
the  pressure  drop  is  significantly  higher  at  the  channels  that 
carry  more  fluid.  As  the  density  is  conserved,  this  drive  is  ex¬ 
pressed  in  higher  velocities  through  those  channels,  as  shown 
in  Fig.  14.  Both  figures  show  a  physical  meaningful  fluid 
field. 
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Fig.  12.  Average  velocities  values  calculated  from  the  experiments,  compared  with  those  obtained  by  numerical  simulation  using  the  138  x  138  grid  and  50,000 
iteration  steps. 
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Fig.  13.  Pressure  map  obtained  from  the  numerical  simulation  with  the  same  experimental  conditions  as  in  Fig.  12. 
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Fig.  14.  The  corresponding  velocity  map  obtained  from  the  numerical  simulation  for  the  conditions  in  Fig.  12. 


64 


F.  Barreras  et  al.  /  Journal  of  Power  Sources  144  (2005)  54-66 


Fig.  15  shows  a  zoom  of  the  region  near  the  entrance, 
showing  velocity  (left)  and  pressure  fields  with  streamlines 
(right).  Streamlines  are  used  to  feature  the  position  of  the 
eddy,  which,  as  mentioned  above,  has  a  large  influence  on 
the  distribution  of  the  flow  in  the  bipolar  plate.  The  back¬ 
ground  depicts  the  pressure  isocontour  levels  and  velocity 
modulus,  respectively.  As  it  has  been  discussed  in  the  exper¬ 
imental  results,  it  can  also  be  observed  the  existence  of  recir¬ 
culation  bubbles  at  the  entrance  of  the  horizontal  channels. 
Although  most  of  the  fluid  tends  to  flow  down  through  the  lat¬ 
eral  channel,  a  significant  amount  of  fluid  is  driven  through 
the  first  two  channels.  The  ultimate  reason  is  probably  the 
higher  pressures  originated  near  the  zones  where  the  fluid 
impinges  against  the  bottom  rib  of  the  second  channel.  These 
higher-pressure  areas  near  the  entrance  of  the  channels  force 


the  fluid  to  flow  preferably  through  those  channels.  Fig.  16 
shows  the  pressure  and  velocity  field,  as  in  Fig.  12,  for  the 
first  channels  but  at  the  zone  opposite  to  the  entrance. 

Figs.  17  and  18  show  the  behavior  of  the  velocity  and 
pressure  at  the  bottom  of  the  bipolar  plate.  It  is  important  to 
note  the  increase  in  the  pressure  drop  of  the  two  last  channels, 
which  can  also  be  verified  analyzing  the  numerical  values 
depicted  in  Fig.  13.  This  fact  motivates  that  a  large  amount 
of  fluid  flows  through  channels  17  and  18  causing  an  increase 
of  the  velocity  in  this  part  of  the  plate.  At  the  same  time,  it 
should  also  be  noted  that  there  is  no  recirculation  bubble  (or 
it  is  very  weak)  in  the  entrance  of  these  channels.  It  is  evident 
that  this  complex  flow  behavior  is  both  studied  and  explained 
in  an  easy  way  using  numerical  simulation  analysis.  A  slight 
shifting  of  the  conditions  in  this  area  may  cause  an  increase 


|V|  (m/s) 

|  1.3540E+00 
2.201 8E-01 

1  8.1 428E-02 
4.1 937E-02 

2.481 2E-02 

I 

3.2040E-06 


P  (Pa) 

3.2461  E+03 
2.5783E+03 
2.3796E+03 
2.2877E+03 
2.1842E+03 
1.3118E+03 


Fig.  15.  Zoom  of  the  region  near  the  entrance  of  the  flow  at  the  bipolar  plate  showing  the  velocity  (a)  and  pressure  fields  (b).  Streamlines  depict  the  position 
of  the  eddies. 


|V|  (m/s) 
1.3540E+00 

2.2018E-01 

I  8.1428E-02 
4.1937E-02 
2.481 2E-02 
3.2040E-06 


P  (Pa) 

2.691 4E+03 
2.4445E+03 
2.3139E+03 
2.1493E+03 
2.1412E+03 
2.1 343 E+03 


Fig.  16.  Velocity  (a)  and  pressure  (b)  fields  for  the  first  channels  to  the  zone  corresponding  to  the  discharge  of  flow  to  the  left  lateral  channel. 
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Fig.  17.  Zoom  image  showing  the  velocity  (a)  and  pressure  (b)  fields  for  the  bottom  channels  to  the  zone  corresponding  to  the  entrance  of  flow. 
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Fig.  18.  Zoom  of  the  region  close  to  the  exit  of  the  flow  to  the  bipolar  plate  showing  the  velocity  (a)  and  pressure  fields  (b). 


of  the  amount  of  flow  through  this  channel  and,  as  mentioned 
in  the  Section  5,  may  be  the  cause  for  the  slower  convergence 
of  the  velocity  in  that  channel. 

It  should  be  pointed  out  that  the  deficient  two-dimensional 
flow  distribution  observed  in  both  experiments  and  numeri¬ 
cal  simulations  is  only  representative  of  the  particular  bipolar 
plate  topology  tested  in  this  work.  Even  when  similar  non- 
homogeneous  flow  patterns  have  been  previously  verified,  for 
example,  for  the  square  spot  gas  flow  topology  [11,12],  the 
present  results  cannot  be  extrapolated  to  other  plate  geome¬ 
tries.  It  is  also  known  that  with  serpentine  and  interdigitated 
designs,  as  well  as  with  flat  porous  plates;  the  gas  flow  can 
be  more  uniformly  distributed  over  the  whole  electrode  area. 
However,  a  drawback  when  these  bipolar  plates  are  used  is 
that  higher  pressures  are  required  for  the  gas  to  reach  the  elec¬ 


trodes,  and  some  designs  are  difficult  to  manufacture.  These 
disadvantages  are  overcome  by  the  series-parallel  distribu¬ 
tion  scheme,  where  a  preferred  flow  path  is  maintained  but 
with  reduced  pressure  drops.  Studies  similar  to  the  one  here 
presented  could  be  performed  with  alternative  plate  models, 
assessing  and  comparing  the  efficiency  of  the  gas  flow  distri¬ 
bution  in  the  cell  plane. 

7.  Conclusions 

Based  on  flow  visualization  and  velocity  measurement  re¬ 
sults,  a  worrying  difference  in  the  velocity  values  inside  a 
parallel  channel  bipolar  plate  has  been  found.  At  the  same 
time,  the  similitude  between  the  dye  patterns  obtained  circu- 
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lating  water  and  a  glycerin-water  mixture  for  the  two  flow 
conditions  tested  demonstrates  a  weak  dependence  on  both 
flow  velocity  and  physical  properties  of  the  fluids.  It  is  evident 
that  the  bipolar  plate  tested  does  not  satisfy  the  requirements 
of  homogeneous  distribution  of  the  flow  over  the  electrodes. 
The  presence  of  channels  along  which  the  fluid  circulates 
very  slowly,  compared  to  the  preferred  paths  through  the  lat¬ 
eral  channels  will  probably  lead  to  an  unbalanced  use  of  the 
catalyst,  and  an  overall  efficiency  lower  than  expected.  This 
non-homogeneity  can  be  partly  explained  by  the  formation 
of  recirculation  bubbles  at  the  channel  entrance  in  the  central 
zone  of  the  plate.  This  effect  could  be  alleviated  by  modi¬ 
fying  the  channel  entrance,  avoiding  the  connection  at  right 
angle  with  the  channels.  Simulations  performed  with  a  fi¬ 
nite  volume  code  have  shown  an  excellent  agreement  with 
the  measurements.  In  the  future,  this  tool  could  be  advanta¬ 
geously  used  to  validate  design  modifications. 

At  least  in  the  tested  conditions,  it  seems  evident  that  the 
design  of  this  particular  bipolar  plate  could  be  substantially 
improved,  optimizing  the  distribution  of  the  reactant  gases 
inside  the  fuel  cell. 

Acknowledgement 

This  research  has  been  partially  supported  by  the  Network 
of  Fuel  Cell  Batteries  of  the  Spanish  Council  for  Scientific 
Research. 

References 

[1]  Hydrogen,  Energy  and  Fuel  Cells — A  vision  for  our  future,  High 
Level  Group  for  Hydrogen  and  Fuel  Cells,  Summary  Report,  June 
2003. 

[2]  V.  Gurau,  H.  Liu,  S.  Kakag,  Two-dimensional  model  for  proton  ex¬ 
change  membrane  fuel  cells,  AIChE  J.  44  (11)  (1998)  2410-2422. 

[3]  R  Costamagna,  Transport  phenomena  in  polymeric  membrane  fuel 
cells,  Chem.  Eng.  Sci.  56  (2001)  323-332. 


[4]  L.  Wang,  A.  Husar,  T.  Zhou,  H.  Liu,  A  parametric  study  of  PEMFC 
performance,  Int.  J.  Hydrogen  Energy  28  (2003)  1263-1272. 

[5]  A.B.  Geiger,  A.  Tsukada,  E.  Lehmann,  P.  Vontobel,  G.G.  Scherer,  In 
situ  investigation  of  two-phase  flow  patterns  in  flow  fields  of  PEFC’s 
using  neutron  radiography,  Fuel  Cells  2  (2)  (2002)  92-98. 

[6]  G.  Bender,  M.S.  Wilson,  T.A.  Zawodzinski,  Further  refinements  in 
the  segmented  cell  approach  to  diagnosing  performance  in  polymer 
electrolyte  fuel  cells,  J.  Power  Sources  123  (2)  (2003)  163-171. 

[7]  K.  Tuber,  D.  Pocza,  Ch.  Hebling,  Visualization  of  water  buildup  in 
the  cathode  of  a  transparent  PEM  fuel  cell,  J.  Power  Sources  124 
(2)  (2003)  403-414. 

[8]  M.M.  Mench,  Q.L.  Dong,  C.Y.  Wang,  In  situ  water  distribution  mea¬ 
surements  in  a  polymer  electrolyte  fuel  cell,  J.  Power  Sources  124 
(1)  (2003)  90-98. 

[9]  T.  Bewer,  T.  Beckman,  H.  Dohle,  J.  Mergel,  D.  Stolten,  Novel 
method  for  investigation  of  two-phase  flow  in  liquid  feed  direct 
methanol  fuel  cells  using  an  aqueous  H2O2  solution,  J.  Power 
Sources  125  (1)  (2004)  1-9. 

[10]  V.  Mehta,  J.S.  Cooper,  Review  and  analysis  of  PEM  fuel  cell  design 
and  manufacturing,  J  Power  Sources  114  (2003)  32-53. 

[11]  L.  Carrette,  K.A.  Friedrich,  U.  Stimming,  Fuel  Cells-fundamentals 
and  applications,  Fuel  Cells  1  (1)  (2001)  5-39. 

[12]  P.  Costamagna,  S.  Srinivasan,  Quantum  jumps  in  the  PEMFC  science 
and  technology  from  the  1960s  to  the  year.  Part  II.  Engineering, 
technology  development  and  application  aspects,  J.  Power  Sources 
102  (2001)  253-269. 

[13]  E.  Hontanon,  M.J.  Escudero,  C.  Bautista,  P.  Garcia- Ybarra,  L.  Daza, 
Optimization  of  flow-field  in  polymer  electrolyte  membrane  fuel  cells 
using  computational  fluid  dynamics  techniques,  J.  Power  Sources  86 
(1-2)  (2000)  363-368. 

[14]  H.  Dohle,  R.  Jung,  N.  Kimiaie,  J.  Mergel,  M.  Muller,  Interaction 
between  the  difussion  layer  and  the  flow  field  of  polymer  electrolite 
fue  cells  -  experiments  and  numerical  simulations,  J.  Power  Sources 
124  (2)  (2003)  371-384. 

[15]  Personal  communication  with  researchers  of  the  National  Institute 
of  Aerospatial  Technique  (I.N.T.A.)  of  Spain. 

[16]  R.H.  Perry,  D.W.  Green,  J.O.  Maloney,  Perry’s  Chemical  Engineers 
Handbook,  Sixth  ed.,  McGraw-Hill,  1984. 

[17]  F.M.  White,  Fluid  Mechanics,  McGraw  Hill,  USA,  1979. 

[18]  U.  Brackmann,  Lambdachrome  Laser  Dyes,  Lambda  Physik  GmbH, 
Gottingen,  Germany,  1994,  pp.  160-161. 

[19]  J.H.  Ferziger,  M.  Peric,  Computational  Methods  for  Fluid  Dynamics. 
Springer- Verlag,  2002. 

[20]  S.V.  Patankar,  Numerical  Heat  Transfer  and  Fluid  Flow,  McGraw- 
Hill,  New  York,  USA,  1980. 


